# Port of error_cov.f90 and correl.f90

module PriorCovariance

using LinearAlgebra
using Dates
using Printf
using Statistics
using ..CoreTypes: Domain
using ..Settings: Config
using ..Tools: grid_area
using ..ForwardModel: aggregate_flux_to_domain, PriorFluxes, FluxField, nearest_time_index
using ..TransportInit: initialize_dates, DateInfo

import ..NUMERICAL_SCALE
export CovarianceMatrix, build_prior_covariance, spatial_correlation, temporal_correlation

const NUMSCALE_AREA_TOL = 1.0e-12

function build_numscale_area_scaling(nbox::Int, nt_total::Int;
                                     area::Union{Vector{Float64},Nothing}=nothing,
                                     override::Union{Vector{Float64},Nothing}=nothing)
    nt_total = max(nt_total, 1)
    physical_len = nbox * nt_total

    if override !== nothing
        if length(override) == physical_len && all(isfinite, override) && all(abs.(override) .> NUMSCALE_AREA_TOL)
            scaling = copy(override)
            inv_scaling = 1.0 ./ scaling
            has_scaling = any(abs.(scaling .- 1.0) .> NUMSCALE_AREA_TOL)
            return scaling, inv_scaling, has_scaling
        end
    end

    if area !== nothing && length(area) == nbox && all(isfinite, area) && all(area .> 0)
        base = area ./ NUMERICAL_SCALE
        scaling = repeat(base, nt_total)
        inv_scaling = 1.0 ./ scaling
        has_scaling = any(abs.(scaling .- 1.0) .> NUMSCALE_AREA_TOL)
        return scaling, inv_scaling, has_scaling
    end

    scaling = ones(Float64, physical_len)
    return scaling, scaling, false
end

function search_fortran_dirs(base_dir::AbstractString)
    clean_base = rstrip(base_dir, ['/', '\\'])
    dirs = String[]
    if !isempty(clean_base) && isdir(clean_base)
        push!(dirs, clean_base)
    end
    parent = dirname(clean_base)
    if parent != clean_base
        candidate = joinpath(parent, "FIVOUT_05deg")
        if isdir(candidate) && candidate ∉ dirs
            push!(dirs, candidate)
        end
        candidate = joinpath(parent, "FIVOUT_05deg_debug")
        if isdir(candidate) && candidate ∉ dirs
            push!(dirs, candidate)
        end
        if isdir(parent)
            prev_dirs = String[]
            for entry in readdir(parent; join=true)
                name = basename(entry)
                if startswith(name, "FIVOUT_05deg_prev") && isdir(entry) && entry ∉ dirs
                    push!(prev_dirs, entry)
                end
            end
            sort!(prev_dirs)
            append!(dirs, prev_dirs)
        end
    end
    return dirs
end

function locate_fortran_file(base_dir::AbstractString, filename::AbstractString)
    dirs = search_fortran_dirs(base_dir)
    prefer_prev = get(ENV, "FLEXINVERT_USE_FORTRAN_H", "0") in ("1", "true", "TRUE")
    if prefer_prev
        preferred = [d for d in dirs if occursin("preRepoRun", basename(d))]
        remaining_prev = [d for d in dirs if occursin("FIVOUT_05deg_prev", basename(d)) && !occursin("preRepoRun", basename(d))]
        sort!(remaining_prev; rev=true)
        other_dirs = [d for d in dirs if !occursin("FIVOUT_05deg_prev", basename(d))]
        dirs = vcat(preferred, remaining_prev, other_dirs)
    end
    seen = Set{String}()
    for dir in dirs
        path = joinpath(dir, filename)
        if isfile(path) && !(path in seen)
            push!(seen, path)
            return path
        end
    end
    return nothing
end

"""
    CovarianceMatrix

Represents a prior error covariance matrix with eigenvalue decomposition for efficient operations.
"""
struct CovarianceMatrix
    # Full covariance matrix (stored if small enough)
    full_matrix::Union{Matrix{Float64}, Nothing}

    # Eigenvalue decomposition for efficient operations
    eigenvalues::Vector{Float64}
    eigenvectors::Matrix{Float64}

    # Temporal correlation matrices (from Fortran's error_cov.f90)
    temporal_corr::Matrix{Float64}
    temporal_corr_inv::Matrix{Float64}
    temporal_corr_sqrt::Matrix{Float64}
    temporal_corr_sqrt_inv::Matrix{Float64}

    # Truncation info
    n_modes::Int
    total_variance::Float64
    explained_variance::Float64

    # Spatial dimensions
    nx::Int
    ny::Int
    nt::Int
    nbox::Int

    # Optional physical-space scaling (legacy physical scaling)
    row_scaling::Vector{Float64}
    inv_row_scaling::Vector{Float64}
    has_row_scaling::Bool
    # NUMSCALE·area scaling (primary physical scaling used in Fortran)
    numscale_area_scaling::Vector{Float64}
    inv_numscale_area_scaling::Vector{Float64}
    has_numscale_area_scaling::Bool
    physical_scaling::Vector{Float64}
    inv_physical_scaling::Vector{Float64}
    has_physical_scaling::Bool
    hbht_row_scaling::Vector{Float64}
    inv_hbht_row_scaling::Vector{Float64}
    has_hbht_row_scaling::Bool
    spatial_covariance::Union{Matrix{Float64},Nothing}
    ntstate::Int
    ndt::Int
    fortran_correlation::Union{Matrix{Float64},Nothing}
    fortran_xerr::Union{Vector{Float64},Nothing}
    fortran_area::Union{Vector{Float64},Nothing}
    fortran_err_box::Union{Array{Float64,3},Nothing}
    fortran_spatial_evals::Union{Vector{Float64},Nothing}
    fortran_spatial_evecs::Union{Matrix{Float64},Nothing}
    globerr::Float64
    has_embedded_temporal_basis::Bool
end

"""
    build_prior_covariance(cfg::Config, domain::Domain) -> CovarianceMatrix

Construct prior covariance matrix with spatial and temporal correlations.

# Arguments
- `cfg`: Configuration settings with correlation parameters
- `domain`: Domain information with grid structure

# Returns
- CovarianceMatrix with eigenvalue decomposition for efficient operations
"""
function build_prior_covariance(cfg::Config, domain::Domain)
    return build_prior_covariance(cfg, domain, nothing)
end

function build_prior_covariance(cfg::Config, domain::Domain, fluxes::Union{PriorFluxes,Nothing})
    @info "Building prior covariance matrix"

    nx = domain.nxregrid
    ny = domain.nyregrid
    nbox = domain.nbox

    if !isempty(strip(cfg.fortran_output_path))
        cov_from_fortran = try_load_fortran_covariance(cfg, domain)
        if cov_from_fortran !== nothing
            @info "Loaded prior covariance from Fortran reference artefacts" path=cfg.fortran_output_path
            return cov_from_fortran
        else
            @warn "Failed to import prior covariance from Fortran artefacts; falling back to Julia construction" path=cfg.fortran_output_path
        end
    end

    date_info = initialize_dates(cfg)
    ntstate = max(date_info.ntstate, 1)
    ndt = max(date_info.ndt, 1)
    nt_total = ntstate * ndt

    ndvar = nbox * nt_total

    @info "Covariance dimensions: $nbox spatial boxes, $nt_total time steps = $ndvar total variables"

    # Build spatial correlation matrix
    spatial_corr = spatial_correlation(cfg, domain)

    # Build temporal correlation matrix
    temporal_corr_state = temporal_correlation(cfg, ntstate; stat_step=cfg.statres)
    temporal_corr_subdaily = ndt == 1 ? Matrix{Float64}(I, 1, 1) : Matrix{Float64}(I, ndt, ndt)
    temporal_corr = kron(temporal_corr_state, temporal_corr_subdaily)

    # Build temporal correlation matrices for transformations (matching Fortran's error_cov.f90)
    temporal_state_corr = build_temporal_correlation_matrices(cfg, ntstate; stat_step=cfg.statres)

    # Build full covariance matrix
    error_variances_base = get_error_variances(cfg, domain, fluxes, date_info)
    error_variances = repeat(error_variances_base, ntstate)

    cov_matrix = build_full_covariance(cfg, domain, spatial_corr, temporal_corr, error_variances)

    # Perform eigenvalue decomposition
    @info "Performing eigenvalue decomposition"
    eigenvals, eigenvecs = eigen(Symmetric(cov_matrix))

    # Sort eigenvalues in descending order
    perm = sortperm(eigenvals, rev=true)
    eigenvals = eigenvals[perm]
    eigenvecs = eigenvecs[:, perm]

    # Calculate variance statistics
    total_var = sum(eigenvals)

    # Determine number of modes to retain
    n_modes = determine_n_modes(eigenvals, cfg)

    # Truncate to significant modes
    eigenvals_trunc = eigenvals[1:n_modes]
    eigenvecs_trunc = eigenvecs[:, 1:n_modes]

    # Apply ridge regularisation to avoid near-singular modes
    ridge_fraction = 1e-6
    largest_eigen = maximum(eigenvals_trunc)
    ridge_floor = max(largest_eigen * ridge_fraction, 1e-12)
    eigenvals_trunc = max.(eigenvals_trunc, ridge_floor)

    explained_var = sum(eigenvals_trunc) / total_var

    @info "Retained $n_modes modes explaining $(round(explained_var*100, digits=1))% of variance"

    # Store full matrix only if small
    full_matrix = nothing
    if ndvar <= 1000  # Threshold for storing full matrix
        full_matrix = cov_matrix
    end

    physical_len = ndvar
    numscale_area_scaling, inv_numscale_area_scaling, has_numscale_area_scaling =
        build_numscale_area_scaling(nbox, nt_total; area=domain.area_box)
    physical_scaling = ones(Float64, physical_len)
    inv_physical_scaling = ones(Float64, physical_len)
    has_physical_scaling = false
    row_scaling = ones(Float64, physical_len)
    inv_row_scaling = ones(Float64, physical_len)
    has_row_scaling = false
    hbht_row_scaling = ones(Float64, physical_len)
    inv_hbht_row_scaling = ones(Float64, physical_len)
    has_hbht_row_scaling = false

    return CovarianceMatrix(
        full_matrix,
        eigenvals_trunc,
        eigenvecs_trunc,
        temporal_state_corr.cort,
        temporal_state_corr.icort,
        temporal_state_corr.sqcort,
        temporal_state_corr.isqcort,
        n_modes,
        total_var,
        explained_var,
        nx, ny, nt_total, nbox,
        row_scaling,
        inv_row_scaling,
        has_row_scaling,
        numscale_area_scaling,
        inv_numscale_area_scaling,
        has_numscale_area_scaling,
        physical_scaling,
        inv_physical_scaling,
        has_physical_scaling,
        hbht_row_scaling,
        inv_hbht_row_scaling,
        has_hbht_row_scaling,
        nothing,
        ntstate,
        ndt,
        nothing,
        nothing,
        nothing,
        nothing,
        nothing,
        nothing,
        cfg.globerr,
        true
    )
end

function read_fortran_lines(path::AbstractString)
    rows = Vector{Vector{Float64}}()
    open(path, "r") do io
        for raw in eachline(io)
            line = strip(raw)
            isempty(line) && continue
            tokens = split(line)
            parsed = Float64[]
            for token in tokens
                sanitized = replace(token, 'D' => 'E', 'd' => 'E')
                push!(parsed, parse(Float64, sanitized))
            end
            push!(rows, parsed)
        end
    end
    return rows
end

function read_fortran_vector(path::AbstractString)
    rows = read_fortran_lines(path)
    isempty(rows) && return Float64[]
    lengths = unique(length.(rows))
    if length(lengths) != 1
        error("Inconsistent column counts in vector file $path")
    end
    ncols = first(lengths)
    if ncols == 1
        # Single column: return as-is
        return [row[1] for row in rows]
    elseif ncols == 2 && all(abs(row[1] - i) <= 1e-6 for (i, row) in enumerate(rows))
        # Two columns with sequential indices: extract values only
        debug_enabled = get(ENV, "FLEXINVERT_DEBUG_OUTPUT", "0") in ("1", "true", "TRUE")
        if debug_enabled
            @info "Detected (index, value) format in Fortran vector" path n_entries=length(rows)
        end
        return [row[2] for row in rows]
    elseif length(rows) == 1
        # Single row with multiple values
        return copy(rows[1])
    else
        # Multiple rows, flatten all values
        flattened = Float64[]
        for row in rows
            append!(flattened, row)
        end
        return flattened
    end
end

function read_fortran_matrix(path::AbstractString)
    rows = read_fortran_lines(path)
    isempty(rows) && return zeros(Float64, 0, 0)
    ncols = length(first(rows))
    if any(length(row) != ncols for row in rows)
        error("Inconsistent column counts in matrix file $path")
    end
    mat = Matrix{Float64}(undef, length(rows), ncols)
    for (i, row) in enumerate(rows)
        for (j, val) in enumerate(row)
            mat[i, j] = val
        end
    end
    return mat
end

function load_area_box_fortran(domain::Domain, fortran_dir::AbstractString, nbox::Int)
    if length(domain.area_box) == nbox && all(isfinite, domain.area_box)
        return copy(domain.area_box)
    end

    area_path = locate_fortran_file(fortran_dir, "area_box.txt")
    if area_path !== nothing && isfile(area_path)
        area_vec = read_fortran_vector(area_path)
        if length(area_vec) == nbox
            return area_vec
        else
            @warn "area_box.txt length mismatch" path=area_path expected=nbox actual=length(area_vec)
        end
    else
        @warn "area_box.txt not found; falling back to domain metadata" path=area_path
    end

    if length(domain.area_box) == nbox
        return copy(domain.area_box)
    end

    @warn "Domain area_box length mismatch; using zeros" expected=nbox actual=length(domain.area_box)
    return zeros(Float64, nbox)
end

function build_spatial_covariance_from_fortran(cfg::Config, corr::Matrix{Float64},
                                               xerr::Vector{Float64}, area_box::Vector{Float64},
                                               ndt::Int)
    nbox = size(corr, 1)
    ndvar = nbox * max(ndt, 1)
    length(xerr) == ndvar || error("xerr length $(length(xerr)) incompatible with nbox=$nbox, ndt=$ndt")

    cov = zeros(Float64, ndvar, ndvar)
    covsum = 0.0
    inv_ndt = ndt > 0 ? Float64(ndt) : 1.0

    for nt = 1:max(ndt, 1)
        row_offset = (nt - 1) * nbox
        for j = 1:nbox
            area_j = j <= length(area_box) ? area_box[j] : 0.0
            idx_j = row_offset + j
            xerr_j = xerr[idx_j]
            for i = 1:nbox
                idx_i = row_offset + i
                val = corr[i, j] * xerr[idx_i] * xerr_j
                cov[idx_i, idx_j] = val
                area_i = i <= length(area_box) ? area_box[i] : 0.0
                covsum += val * area_i * area_j
            end
        end
    end

    covsum /= inv_ndt^2
    scale_factor = 1.0
    if !cfg.lognormal && cfg.globerr > 0 && covsum > 0.0
        toterr = sqrt(covsum) * 3600.0 * 24.0 * 365.0 / 1e9 / NUMERICAL_SCALE
        if toterr > 0.0 && isfinite(cfg.globerr)
            scale_factor = cfg.globerr / toterr
            cov .*= scale_factor^2
        end
    end

    return Symmetric(cov), scale_factor
end

function try_load_fortran_covariance(cfg::Config, domain::Domain)
    fortran_dir = strip(cfg.fortran_output_path)
    isempty(fortran_dir) && return nothing
    isdir(fortran_dir) || return nothing

    evals_path = joinpath(fortran_dir, "evals.txt")
    evecs_path = joinpath(fortran_dir, "evecs.txt")
    cort_path = joinpath(fortran_dir, "cort.txt")
    sqcort_path = joinpath(fortran_dir, "sqcort.txt")
    icort_path = joinpath(fortran_dir, "icort.txt")

    if !isfile(evals_path) || !isfile(evecs_path)
        return nothing
    end

    eigenvalues = read_fortran_vector(evals_path)
    eigenvectors = read_fortran_matrix(evecs_path)

    date_info = initialize_dates(cfg)
    expected_ntstate = max(date_info.ntstate, 1)
    expected_ndt = max(date_info.ndt, 1)

    temporal_corr_state = isfile(cort_path) ? read_fortran_matrix(cort_path) : Matrix{Float64}(I, expected_ntstate, expected_ntstate)
    if size(temporal_corr_state, 1) == 0
        temporal_corr_state = Matrix{Float64}(I, expected_ntstate, expected_ntstate)
    end
    sqcort_state = isfile(sqcort_path) ? read_fortran_matrix(sqcort_path) : copy(temporal_corr_state)
    icort_state = if isfile(icort_path)
        read_fortran_matrix(icort_path)
    else
        Array{Float64,2}(inv(Symmetric(temporal_corr_state)))
    end

    ntstate = size(temporal_corr_state, 1)
    if ntstate != expected_ntstate
        @warn "Temporal state count from Fortran artefact differs from config" artefact_ntstate=ntstate config_ntstate=expected_ntstate
    end
    ntstate = max(ntstate, 1)

    corr_path = locate_fortran_file(fortran_dir, "correl.txt")
    xerr_path = locate_fortran_file(fortran_dir, "xerr.txt")
    err_box_path = locate_fortran_file(fortran_dir, "err_box.txt")

    spatial_cov = nothing
    ndt = expected_ndt
    nbox_cov = domain.nbox
    corr_matrix = nothing
    xerr_values = nothing
    area_box_vals = length(domain.area_box) == domain.nbox ? copy(domain.area_box) : nothing
    err_box_array = nothing

    if corr_path !== nothing && xerr_path !== nothing
        corr_matrix = read_fortran_matrix(corr_path)
        nbox_cov = size(corr_matrix, 1)
        if size(corr_matrix, 2) != nbox_cov
            @warn "correl.txt is not square; falling back to eigen reconstruction" path=corr_path
        else
            xerr_values = read_fortran_vector(xerr_path)
            if isempty(xerr_values)
                @warn "xerr.txt empty; falling back to eigen reconstruction" path=xerr_path
            else
                ndt_from_xerr = max(1, Int(round(length(xerr_values) / max(nbox_cov, 1))))
                if nbox_cov * ndt_from_xerr != length(xerr_values)
                    @warn "xerr length mismatch; expected multiples of nbox" length=length(xerr_values) nbox=nbox_cov
                else
                    area_box_vals = load_area_box_fortran(domain, fortran_dir, nbox_cov)
                    spatial_cov_sym, _ = build_spatial_covariance_from_fortran(cfg, corr_matrix, xerr_values, area_box_vals, ndt_from_xerr)
                    spatial_cov = Matrix(spatial_cov_sym)
                    ndt = ndt_from_xerr
                end
            end
        end
    end
    if err_box_path !== nothing && isfile(err_box_path)
        err_matrix = read_fortran_matrix(err_box_path)
        reshaped = reshape_err_box_matrix(err_matrix, nbox_cov, ntstate)
        if reshaped === nothing
            @warn "Failed to reshape err_box.txt; ignoring file" path=err_box_path
        else
            err_box_array, ndt_from_err = reshaped
            ndt = ndt_from_err
        end
    end

    if spatial_cov === nothing
        n_spatial = size(eigenvectors, 1)
        if n_spatial == 0
            @warn "Fortran eigenvector file empty"
            return nothing
        end
        if domain.nbox != 0 && domain.nbox != n_spatial
            @warn "Spatial dimension mismatch between Fortran artefact and domain" fortran_nbox=n_spatial domain_nbox=domain.nbox
        end
        perm = sortperm(eigenvalues, rev=true)
        eigenvalues = eigenvalues[perm]
        eigenvectors = eigenvectors[:, perm]
        spatial_cov = eigenvectors * Diagonal(eigenvalues) * eigenvectors'
        nbox_cov = n_spatial
        ndt = expected_ndt
    end

    ndvar = size(spatial_cov, 1)
    nt_total = ntstate * ndt
    physical_len = ndvar * ntstate

    temporal_corr = Matrix{Float64}(temporal_corr_state)
    temporal_corr_inv = Matrix{Float64}(icort_state)
    temporal_corr_sqrt = Matrix{Float64}(sqcort_state)
    temporal_corr_sqrt_inv = try
        Matrix(inv(Symmetric(sqcort_state)))
    catch err
        @warn "Failed to invert sqcort matrix from Fortran artefacts; recomputing via eigen decomposition" error=err
        temporal_eig = eigen(Symmetric(temporal_corr_state))
        temporal_eig.vectors * Diagonal(1.0 ./ sqrt.(temporal_eig.values)) * temporal_eig.vectors'
    end

    cov_full_sym = Symmetric(kron(temporal_corr_state, spatial_cov))
    eig_full = eigen(cov_full_sym)
    eigenvals_full = eig_full.values
    eigenvecs_full = eig_full.vectors
    eigenvals_full = map(x -> (x < 0 && abs(x) < 1e-12) ? 0.0 : x, eigenvals_full)
    perm_full = sortperm(eigenvals_full, rev=true)
    eigenvals_full = eigenvals_full[perm_full]
    eigenvecs_full = eigenvecs_full[:, perm_full]

    total_var = sum(eigenvals_full)
    n_modes = determine_n_modes(eigenvals_full, cfg)
    eigenvals_trunc = eigenvals_full[1:n_modes]
    eigenvecs_trunc = eigenvecs_full[:, 1:n_modes]
    explained_var = total_var > 0 ? sum(eigenvals_trunc) / total_var : 1.0

    physical_scaling_vec, physical_override_present = load_fortran_physical_scaling(fortran_dir, physical_len)
    numscale_area_scaling, inv_numscale_area_scaling, has_numscale_area_scaling =
        build_numscale_area_scaling(nbox_cov, ntstate * ndt;
                                    area=area_box_vals,
                                    override=physical_override_present ? physical_scaling_vec : nothing)
    if physical_override_present
        physical_scaling_vec = ones(Float64, physical_len)
        inv_physical_scaling_vec = ones(Float64, physical_len)
        has_physical_scaling = false
    else
        inv_physical_scaling_vec = 1.0 ./ physical_scaling_vec
        has_physical_scaling = any(abs.(physical_scaling_vec .- 1.0) .> NUMSCALE_AREA_TOL)
    end
    row_scaling_vec, inv_row_scaling_vec, has_row_scaling = load_fortran_row_scaling(fortran_dir, physical_len)
    hbht_row_scaling_vec, inv_hbht_row_scaling_vec, has_hbht_row_scaling = load_fortran_hbht_scaling(fortran_dir, physical_len)

    return CovarianceMatrix(
        Matrix(cov_full_sym),
        eigenvals_trunc,
        eigenvecs_trunc,
        temporal_corr,
        temporal_corr_inv,
        temporal_corr_sqrt,
        temporal_corr_sqrt_inv,
        n_modes,
        total_var,
        explained_var,
        domain.nxregrid,
        domain.nyregrid,
        nt_total,
        nbox_cov,
        row_scaling_vec,
        inv_row_scaling_vec,
        has_row_scaling,
        numscale_area_scaling,
        inv_numscale_area_scaling,
        has_numscale_area_scaling,
        physical_scaling_vec,
        inv_physical_scaling_vec,
        has_physical_scaling,
        hbht_row_scaling_vec,
        inv_hbht_row_scaling_vec,
        has_hbht_row_scaling,
        spatial_cov === nothing ? nothing : Matrix(spatial_cov),
        size(temporal_corr_state, 1),
        ndt,
        nothing,
        nothing,
        area_box_vals,
        err_box_array,
        length(eigenvalues) > 0 ? copy(eigenvalues) : nothing,
        size(eigenvectors, 1) > 0 ? Matrix(eigenvectors) : nothing,
        cfg.globerr,
        false
    )
end


function reshape_err_box_matrix(data::Matrix{Float64}, nbox::Int, ntstate::Int)
    total_rows = size(data, 1)
    nbox == size(data, 2) || return nothing
    ntstate = max(ntstate, 1)
    total_rows % ntstate == 0 || return nothing
    ndt = Int(total_rows ÷ ntstate)
    arr = Array{Float64,3}(undef, nbox, ntstate, ndt)
    row = 1
    @inbounds for nts in 1:ntstate
        for bin in 1:ndt
            arr[:, nts, bin] .= data[row, :]
            row += 1
        end
    end
    return arr, ndt
end

function load_fortran_physical_scaling(dir::AbstractString, expected_len::Int)
    base_dir = rstrip(dir, [Char(0x2f), Char(0x5c)])
    parent_dir = dirname(base_dir)
    search_dirs = String[base_dir]
    debug_dir = joinpath(parent_dir, "FIVOUT_05deg_debug")
    if isdir(debug_dir)
        push!(search_dirs, debug_dir)
    end

    debug_enabled = get(ENV, "FLEXINVERT_DEBUG_OUTPUT", "0") in ("1", "true", "TRUE")

    for base in search_dirs
        exact = joinpath(base, "physical_scaling_vector.txt")
        if isfile(exact)
            vec = read_fortran_vector(exact)
            if length(vec) == expected_len && all(isfinite, vec) && all(abs.(vec) .> 0.0)
                max_abs = maximum(abs.(vec))
                if max_abs <= 1.0e9
                    @info "Loaded physical scaling vector" path=exact n_entries=expected_len min=minimum(abs.(vec)) max=max_abs median=median(vec)
                    # Cache parsed vector for auditing when debug output enabled
                    if debug_enabled
                        output_dir = get(ENV, "FLEXINVERT_OUTPUT_DIR", "")
                        if !isempty(output_dir) && isdir(output_dir)
                            cache_path = joinpath(output_dir, "physical_scaling_vector_parsed.txt")
                            try
                                open(cache_path, "w") do io
                                    for val in vec
                                        @printf(io, "%.16e\n", val)
                                    end
                                end
                                @info "Cached parsed physical scaling vector" cache_path
                            catch e
                                @warn "Failed to cache physical scaling vector" cache_path exception=e
                            end
                        end
                    end
                    return vec, true
                else
                    @warn "Physical scaling vector magnitude unrealistic; ignoring file" path=exact max_abs=max_abs
                end
            elseif length(vec) == expected_len
                @warn "Physical scaling vector contains zero entries; ignoring file" path=exact
            else
                @warn "Physical scaling vector length mismatch on first read" path=exact expected_len actual_len=length(vec)
            end
        end
    end

    candidates = String[]
    for base in search_dirs
        isdir(base) || continue
        for fname in readdir(base; join=false)
            occursin("physical_scaling_vector", fname) || continue
            push!(candidates, joinpath(base, fname))
        end
    end

    if isempty(candidates)
        @info "No physical scaling vector found; using unity scaling" expected_len
        return ones(Float64, expected_len), false
    end
    sort!(candidates)

    for path in candidates
        vec = read_fortran_vector(path)
        if length(vec) == expected_len
            if any(!isfinite, vec)
                @warn "Physical scaling vector contains non-finite values; skipping file" path
                continue
            elseif any(abs.(vec) .== 0.0)
                @warn "Physical scaling vector contains zero entries; skipping file" path
                continue
            end
            max_abs = maximum(abs.(vec))
            if max_abs > 1.0e9
                @warn "Physical scaling vector magnitude unrealistic; skipping file" path max_abs=max_abs
                continue
            end
            @info "Loaded physical scaling vector from fallback" path n_entries=expected_len
            return vec, true
        else
            @warn "Physical scaling vector length mismatch" path expected_len actual_len=length(vec)
        end
    end

    @info "No valid physical scaling vector found; using unity scaling" expected_len
    return ones(Float64, expected_len), false
end
function load_fortran_row_scaling(dir::AbstractString, expected_len::Int)
    base_dir = rstrip(dir, [Char(0x2f), Char(0x5c)])
    parent_dir = dirname(base_dir)
    search_dirs = String[base_dir]
    debug_dir = joinpath(parent_dir, "FIVOUT_05deg_debug")
    if isdir(debug_dir)
        push!(search_dirs, debug_dir)
    end

    for base in search_dirs
        path = joinpath(base, "sym_row_scaling_vector.txt")
        if isfile(path)
            vec = read_fortran_vector(path)
            if length(vec) == expected_len && all(isfinite, vec) && all(abs.(vec) .> 0)
                inv_vec = 1.0 ./ vec
                return vec, inv_vec, true
            end
        end
    end

    ones_vec = ones(Float64, expected_len)
    return ones_vec, ones_vec, false
end

function load_fortran_hbht_scaling(dir::AbstractString, expected_len::Int)
    base_dir = rstrip(dir, [Char(0x2f), Char(0x5c)])
    parent_dir = dirname(base_dir)
    candidates = String[
        joinpath(base_dir, "hbht_scaling_vector.txt"),
        joinpath(parent_dir, "FIVOUT_05deg_debug", "hbht_scaling_vector.txt"),
        joinpath(parent_dir, "FIVOUT_05deg", "hbht_scaling_vector.txt")
    ]

    for path in candidates
        if isfile(path)
            vec = read_fortran_vector(path)
            vals = if size(vec, 2) == 2
                vec[:, 2]
            else
                vec
            end
            if length(vals) == expected_len && all(isfinite, vals) && all(abs.(vals) .> 0)
                inv_vals = 1.0 ./ vals
                return vals, inv_vals, true
            elseif length(vals) != expected_len
                @warn "HBHᵀ scaling vector length mismatch" path expected_len actual_len=length(vals)
            else
                @warn "HBHᵀ scaling vector contains invalid entries; skipping file" path
            end
        end
    end

    ones_vec = ones(Float64, expected_len)
    return ones_vec, ones_vec, false
end

"""
    spatial_correlation(cfg::Config, domain::Domain) -> Matrix{Float64}

Read spatial correlation matrix from correl.txt file (matches Fortran behavior).
If correl.txt is not available, fall back to computing using exponential decay with distance.
"""
function spatial_correlation(cfg::Config, domain::Domain)
    nbox = domain.nbox
    @info "Reading spatial correlation matrix ($nbox x $nbox)"

    # Try to read from correl.txt file first (matches Fortran behavior)
    # Look for correl.txt in both 5° and 0.1° output directories
    # Use absolute path for reliability
    base_dir = "/home/linden/Downloads/FLEXINVERT/run/work/minimal_env/flexinvert_run"

    # Try 5° version first (more likely for current run)
    correl_path = joinpath(base_dir, "FIVOUT_05deg", "correl.txt")
    if !isfile(correl_path)
        @info "5° correl.txt not found at: $correl_path"
        # Fall back to 0.1° version
        correl_path = joinpath(base_dir, "FIVOUT_01deg", "correl.txt")
        @info "Trying 0.1° correl.txt at: $correl_path"
    else
        @info "Found 5° correl.txt at: $correl_path"
    end

    if isfile(correl_path)
        @info "Reading spatial correlation from: $correl_path"
        try
            return read_correl_txt(correl_path, nbox)
        catch err
            @warn "Failed to read correl.txt file; computing spatial correlation from scratch" err
        end
    else
        @info "correl.txt not found at $correl_path; computing spatial correlation from scratch"
    end

    # Fallback: compute spatial correlation matrix using exponential decay with distance
    @info "Computing spatial correlation matrix ($nbox x $nbox)"
    assert_unique_box_mapping(domain)

    corr = zeros(Float64, nbox, nbox)

    # Get correlation length scales
    sigma_land = get_correlation_length(cfg, "land", 300.0)  # km, default
    sigma_ocean = get_correlation_length(cfg, "ocean", 500.0)  # km, default

    # Calculate correlation between all box pairs
    for i in 1:nbox
        # Find grid coordinates for box i
        ix1, jy1, is_land1 = find_box_coordinates(domain, i)
        lon1, lat1 = grid_coordinates(domain, ix1, jy1)
        sigma1 = is_land1 ? sigma_land : sigma_ocean

        for j in i:nbox  # Only upper triangle, use symmetry
            if i == j
                corr[i, j] = 1.0
            else
                # Find grid coordinates for box j
                ix2, jy2, is_land2 = find_box_coordinates(domain, j)
                lon2, lat2 = grid_coordinates(domain, ix2, jy2)
                sigma2 = is_land2 ? sigma_land : sigma_ocean

                # Calculate distance in km
                dist_km = haversine_distance(lat1, lon1, lat2, lon2)

                # Use geometric mean of correlation lengths
                sigma_eff = sqrt(sigma1 * sigma2)

                # Exponential decay correlation
                corr[i, j] = exp(-dist_km / sigma_eff)
                corr[j, i] = corr[i, j]  # Symmetry
            end
        end
    end

    return corr
end

"""
    temporal_correlation(cfg::Config, nt::Int) -> Matrix{Float64}

Calculate temporal correlation matrix (e.g., autoregressive).
"""
function temporal_correlation(cfg::Config, nt::Int; stat_step::Float64=1.0)
    if nt == 1
        return reshape([1.0], 1, 1)
    end

    @info "Computing temporal correlation matrix ($nt x $nt)"

    corr = zeros(Float64, nt, nt)
    step = max(stat_step, 1e-6)

    if hasfield(typeof(cfg), :sigmatime) && getfield(cfg, :sigmatime) > 0
        tau = getfield(cfg, :sigmatime)
        for i in 1:nt
            for j in 1:nt
                lag = abs(i - j)
                corr[i, j] = exp(-lag * step / tau)
            end
        end
    else
        temporal_scale = get_temporal_correlation(cfg, 0.5)
        for i in 1:nt
            for j in 1:nt
                lag = abs(i - j)
                corr[i, j] = temporal_scale^lag
            end
        end
    end

    return corr
end

"""
    build_temporal_correlation_matrices(cfg::Config, nt::Int; stat_step::Float64=1.0) -> NamedTuple

Build temporal correlation matrices matching Fortran's error_cov.f90 implementation.

Returns a NamedTuple with:
- cort: temporal correlation matrix
- icort: inverse of temporal correlation matrix
- sqcort: square root of temporal correlation matrix
- isqcort: inverse of square root of temporal correlation matrix
"""
function build_temporal_correlation_matrices(cfg::Config, nt::Int; stat_step::Float64=1.0)
    if nt == 1
        cort = reshape([1.0], 1, 1)
        icort = copy(cort)
        sqcort = copy(cort)
        isqcort = copy(cort)
        return (cort=cort, icort=icort, sqcort=sqcort, isqcort=isqcort)
    end

    @info "Building temporal correlation matrices ($nt x $nt)"

    # Build time step vector (matches Fortran: tsteps(i) = i * statres)
    tsteps = collect(1.0:1.0:Float64(nt)) .* stat_step

    # Build temporal correlation matrix (matches Fortran logic)
    cort = zeros(Float64, nt, nt)
    for i in 1:nt
        for j in i:nt
            # Calculate time difference (matches Fortran's dtmat)
            dt = abs(tsteps[i] - tsteps[j])

            # Exponential decay correlation (matches Fortran's exp(-dt/sigmatime))
            if hasfield(typeof(cfg), :sigmatime) && cfg.sigmatime > 0
                cort[i, j] = exp(-dt / cfg.sigmatime)
            else
                # Default correlation if sigmatime not specified
                cort[i, j] = 0.5^dt
            end
            cort[j, i] = cort[i, j]  # Symmetry
        end
    end

    # Compute eigenvalue decomposition (matches Fortran's eigen-decomposition)
    eigen_vals, eigen_vecs = eigen(Symmetric(cort))

    # Ensure positive eigenvalues (avoid numerical issues)
    eigen_vals = max.(eigen_vals, 1e-15)

    # Compute inverse of cort: icort = V * diag(1/λ) * V^T
    icort = eigen_vecs * Diagonal(1.0 ./ eigen_vals) * eigen_vecs'

    # Compute square root of cort: sqcort = V * diag(√λ) * V^T
    sqcort = eigen_vecs * Diagonal(sqrt.(eigen_vals)) * eigen_vecs'

    # Compute inverse of square root: isqcort = V * diag(1/√λ) * V^T
    isqcort = eigen_vecs * Diagonal(1.0 ./ sqrt.(eigen_vals)) * eigen_vecs'

    return (cort=cort, icort=icort, sqcort=sqcort, isqcort=isqcort)
end

"""
    build_full_covariance(cfg::Config, domain::Domain, spatial_corr::Matrix{Float64},
                         temporal_corr::Matrix{Float64}) -> Matrix{Float64}

Build full space-time covariance matrix from spatial and temporal components.
"""
function build_full_covariance(cfg::Config, domain::Domain, spatial_corr::Matrix{Float64},
                              temporal_corr::Matrix{Float64}, error_variances::Vector{Float64})
    nbox = domain.nbox
    nt = size(temporal_corr, 1)
    ndvar = nbox * nt

    @info "Building full covariance matrix ($ndvar x $ndvar)"

    cov = zeros(Float64, ndvar, ndvar)

    # Build block structure: temporal correlation within each spatial box pair
    for t1 in 1:nt
        for t2 in 1:nt
            temp_corr = temporal_corr[t1, t2]

            for i in 1:nbox
                for j in 1:nbox
                    spat_corr = spatial_corr[i, j]

                    # Combined correlation
                    combined_corr = spat_corr * temp_corr

                    # Convert to covariance using error standard deviations
                    row_idx = (t1 - 1) * nbox + i
                    col_idx = (t2 - 1) * nbox + j

                    cov[row_idx, col_idx] = combined_corr *
                        sqrt(error_variances[row_idx] * error_variances[col_idx])
                end
            end
        end
    end

    if !cfg.lognormal && cfg.globerr > 0 && !isempty(domain.area_box)
        covsum = 0.0
        nt_float = Float64(nt)
        for t in 1:nt
            base = (t - 1) * nbox
            for j in 1:nbox
                area_j = j <= length(domain.area_box) ? domain.area_box[j] : 0.0
                area_j <= 0 && continue
                for i in 1:nbox
                    area_i = i <= length(domain.area_box) ? domain.area_box[i] : 0.0
                    area_i <= 0 && continue
                    covsum += cov[base + i, base + j] * area_i * area_j
                end
            end
        end
        covsum /= nt_float^2
        if covsum > 0.0
            toterr = sqrt(covsum) * 3600.0 * 24.0 * 365.0 / 1e9 / NUMERICAL_SCALE
            if toterr > 0.0
                erscalar = cfg.globerr / toterr
                if isfinite(erscalar) && erscalar > 0.0
                    cov .*= erscalar^2
                    @info "Scaling prior covariance to match globerr target" total_error_before=toterr scale_factor=erscalar
                end
            end
        end
    end

    return cov
end

# Helper functions

"""
    read_correl_txt(correl_path::String, nbox::Int) -> Matrix{Float64}

Read spatial correlation matrix from Fortran's correl.txt file format.

The correl.txt file contains correlation values in Fortran scientific notation
format with values arranged in rows. Each row contains multiple correlation
values separated by spaces.

# Arguments
- `correl_path`: Path to the correl.txt file
- `nbox`: Expected size of the correlation matrix (nbox x nbox)

# Returns
- Spatial correlation matrix of size (nbox, nbox)
"""
function read_correl_txt(correl_path::String, nbox::Int)
    @info "Reading spatial correlation matrix from correl.txt"

    # Open and read the file
    lines = readlines(correl_path)

    # Parse all floating point values from the file
    values = Float64[]
    for line in lines
        # Skip empty lines
        if isempty(strip(line))
            continue
        end

        # Split line and parse each value
        tokens = split(strip(line))
        for token in tokens
            try
                val = parse(Float64, token)
                push!(values, val)
            catch err
                @warn "Failed to parse token '$token' in correl.txt" err
            end
        end
    end

    # Check we have enough values
    expected_values = nbox * nbox
    if length(values) < expected_values
        error("correl.txt contains $(length(values)) values but expected $expected_values for $nbox x $nbox matrix")
    end

    # Create correlation matrix
    corr = zeros(Float64, nbox, nbox)

    # Fill matrix with parsed values (column-major order like Fortran)
    idx = 1
    for i in 1:nbox
        for j in 1:nbox
            if idx <= length(values)
                corr[i, j] = values[idx]
                idx += 1
            else
                error("Ran out of values while filling correlation matrix at position ($i, $j)")
            end
        end
    end

    # Ensure diagonal is exactly 1.0 (numerical precision)
    for i in 1:nbox
        corr[i, i] = 1.0
    end

    @info "Successfully read $(length(values)) values into $nbox x $nbox correlation matrix"
    return corr
end

"""
    get_correlation_length(cfg::Config, surface_type::String, default::Float64) -> Float64

Get correlation length from configuration.
"""
function get_correlation_length(cfg::Config, surface_type::String, default::Float64)
    if surface_type == "land"
        return hasfield(typeof(cfg), :sigma_land) ? cfg.sigma_land : default
    elseif surface_type == "ocean"
        return hasfield(typeof(cfg), :sigma_ocean) ? cfg.sigma_ocean : default
    else
        return default
    end
end

"""
    get_temporal_correlation(cfg::Config, default::Float64) -> Float64

Get temporal correlation parameter from configuration.
"""
function get_temporal_correlation(cfg::Config, default::Float64)
    return hasfield(typeof(cfg), :temporal_corr) ? cfg.temporal_corr : default
end

"""
    find_box_coordinates(domain::Domain, box_id::Int) -> Tuple{Int, Int, Bool}

Find grid coordinates and land/ocean type for a given box ID.
"""
function find_box_coordinates(domain::Domain, box_id::Int)
    for jy in 1:domain.nyregrid
        for ix in 1:domain.nxregrid
            if domain.nbox_xy[ix, jy] == box_id
                is_land = domain.lsm_box[ix, jy] == 1
                return ix, jy, is_land
            end
        end
    end
    error("Box ID $box_id not found in domain")
end

"""
    assert_unique_box_mapping(domain::Domain)

Ensure each control box is mapped to exactly one underlying grid cell. The legacy
`find_box_coordinates` fast path assumes a bijective grid→box mapping; this guard
prevents silent misuse when redundant assignments exist.
"""
function assert_unique_box_mapping(domain::Domain)
    domain.nbox <= 0 && return
    counts = zeros(Int, domain.nbox)
    for jy in 1:size(domain.nbox_xy, 2)
        for ix in 1:size(domain.nbox_xy, 1)
            box = domain.nbox_xy[ix, jy]
            (1 <= box <= domain.nbox) || continue
            counts[box] += 1
        end
    end
    redundant = findall(>(1), counts)
    isempty(redundant) && return
    details = join([string("box ", b, " has ", counts[b], " cells") for b in redundant], ", ")
    error("Legacy spatial correlation fallback requires unique grid-to-box mapping; " *
          "detected redundant assignments: $details")
end

"""
    grid_coordinates(domain::Domain, ix::Int, jy::Int) -> Tuple{Float64, Float64}

Convert grid indices to longitude/latitude coordinates.
"""
function grid_coordinates(domain::Domain, ix::Int, jy::Int)
    lon = domain.rllx + (ix - 1) * domain.rdx
    lat = domain.rlly + (jy - 1) * domain.rdy
    return lon, lat
end

"""
    haversine_distance(lat1::Float64, lon1::Float64, lat2::Float64, lon2::Float64) -> Float64

Calculate great circle distance between two points in km.
"""
function haversine_distance(lat1::Float64, lon1::Float64, lat2::Float64, lon2::Float64)
    R = 6_371_000.0  # Earth radius in metres

    # Convert to radians
    φ1, φ2 = deg2rad(lat1), deg2rad(lat2)
    Δφ = deg2rad(lat2 - lat1)
    Δλ = deg2rad(lon2 - lon1)

    a = sin(Δφ/2)^2 + cos(φ1) * cos(φ2) * sin(Δλ/2)^2
    c = 2 * atan(sqrt(a), sqrt(1-a))

    return R * c
end

"""
    get_error_variances(cfg::Config, domain::Domain) -> Vector{Float64}

Get error variances for each state variable.
"""
function get_error_variances(cfg::Config, domain::Domain, fluxes::Union{PriorFluxes,Nothing}, date_info::DateInfo)
    nbox = domain.nbox
    nt = max(date_info.ndt, 1)

    if fluxes !== nothing && fluxes.nee !== nothing
        try
            variances = compute_fortran_prior_variances(cfg, domain, fluxes, date_info)
            if variances !== nothing
                return variances
            end
        catch err
            @warn "Fortran-style prior variance construction failed; using heuristic fallback" err
        end
    end

    # Default relative error fallback
    default_rel_error = hasfield(typeof(cfg), :prior_error) ? cfg.prior_error : 0.5

    if fluxes === nothing
        variances = Vector{Float64}(undef, nbox * nt)
        for t in 1:nt
            for box in 1:nbox
                idx = (t - 1) * nbox + box
                area = (box <= length(domain.area_box) && !isempty(domain.area_box)) ? domain.area_box[box] : 1.0
                variances[idx] = (default_rel_error * area)^2
            end
        end
        return variances
    end

    # Heuristic fallback if full port cannot be constructed
    nee_mag = fluxes.nee === nothing ? zeros(Float64, nbox) : mean_abs_flux_per_box(fluxes.nee, domain)
    fff_mag = fluxes.fossil_fuel === nothing ? zeros(Float64, nbox) : mean_abs_flux_per_box(fluxes.fossil_fuel, domain)
    ocn_mag = fluxes.ocean === nothing ? zeros(Float64, nbox) : mean_abs_flux_per_box(fluxes.ocean, domain)

    nee_rel = cfg.flxerr > 0 ? cfg.flxerr : default_rel_error
    fff_rel = cfg.ffferr > 0 ? cfg.ffferr : nee_rel
    ocn_rel = cfg.flxerr > 0 ? cfg.flxerr : default_rel_error

    posterior_std = nee_rel .* nee_mag .+
                    fff_rel .* fff_mag .+
                    (cfg.inc_ocean ? ocn_rel .* ocn_mag : zeros(Float64, nbox))

    minimum_floor = 1e-12
    posterior_std = max.(posterior_std, minimum_floor)

    variances = posterior_std .^ 2
    return repeat(variances, nt)
end

function mean_abs_flux_per_box(flux::FluxField, domain::Domain)
    nbox = domain.nbox
    nx = domain.nxregrid
    ny = domain.nyregrid

    lon_edges = domain_lon_edges(domain)
    lat_edges = domain_lat_edges(domain)

    sums = zeros(Float64, nbox)
    counts = zeros(Int, nbox)

    nt = size(flux.data, 3)
    for t in 1:nt
        slice = flux.data[:, :, t]
        aggregated = aggregate_flux_to_domain(slice, flux.lon, flux.lat, lon_edges, lat_edges)
        aggregated .*= flux.scaling_factor
        for ix in 1:nx
            for iy in 1:ny
                box = domain.nbox_xy[ix, iy]
                if box > 0
                    val = aggregated[ix, iy]
                    if isfinite(val)
                        sums[box] += abs(val)
                        counts[box] += 1
                    end
                end
            end
        end
    end

    for i in 1:nbox
        if counts[i] > 0
            sums[i] /= counts[i]
        else
            sums[i] = 0.0
        end
    end

    return sums
end

function domain_lon_edges(domain::Domain)
    start = domain.rllx
    step = domain.rdx
    nx = domain.nxregrid
    return collect(start:step:(start + step * nx))
end

function domain_lat_edges(domain::Domain)
    start = domain.rlly
    step = domain.rdy
    ny = domain.nyregrid
    return collect(start:step:(start + step * ny))
end

function domain_lon_centers(domain::Domain)
    nx = domain.nxregrid
    centers = Vector{Float64}(undef, nx)
    for ix in 1:nx
        centers[ix] = domain.rllx + (ix - 0.5) * domain.rdx
    end
    return centers
end

function domain_lat_centers(domain::Domain)
    ny = domain.nyregrid
    centers = Vector{Float64}(undef, ny)
    for iy in 1:ny
        centers[iy] = domain.rlly + (iy - 0.5) * domain.rdy
    end
    return centers
end

function precompute_cell_areas(domain::Domain)
    nx = domain.nxregrid
    ny = domain.nyregrid
    lat_centers = domain_lat_centers(domain)
    areas = zeros(Float64, nx, ny)
    for ix in 1:nx
        for iy in 1:ny
            areas[ix, iy] = grid_area(lat_centers[iy], domain.rdy, domain.rdx)
        end
    end
    return areas
end

function area_weighted_box_means(grid_slice::AbstractMatrix{<:Real}, domain::Domain, cell_areas::AbstractMatrix{<:Real})
    nbox = domain.nbox
    sums = zeros(Float64, nbox)
    nx = size(grid_slice, 1)
    ny = size(grid_slice, 2)
    for ix in 1:nx
        for iy in 1:ny
            box = domain.nbox_xy[ix, iy]
            if box > 0
                sums[box] += grid_slice[ix, iy] * cell_areas[ix, iy]
            end
        end
    end
    result = zeros(Float64, nbox)
    for box in 1:nbox
        area = box <= length(domain.area_box) ? domain.area_box[box] : 0.0
        if area > 0
            result[box] = sums[box] / area
        else
            result[box] = 0.0
        end
    end
    return result
end

const HOURS_PER_DAY = 24.0
const SECONDS_PER_HOUR = 3600.0

function write_debug_prior_diagnostics(dir::String, nee_box::Array{Float64,3}, err_box::Array{Float64,3}, variances::Vector{Float64}, flxnee::Array{Float64,3}, nbox::Int, ntstate::Int, nd_nee::Int, ndt::Int)
    mkpath(dir)

    function write_box(path::String, data::Array{Float64,3}, slots::Int)
        open(path, "w") do io
            for state in 1:ntstate
                for slot in 1:slots
                    row = @view data[:, state, slot]
                    for box in 1:nbox
                        @printf(io, "%11.4E%s", row[box], box == nbox ? "\n" : " ")
                    end
                end
            end
        end
    end

    write_box(joinpath(dir, "nee_box.txt"), nee_box, nd_nee)
    write_box(joinpath(dir, "err_box.txt"), err_box, ndt)

    open(joinpath(dir, "xerr.txt"), "w") do io
        for bin in 1:ndt
            offset = (bin - 1) * nbox
            for box in 1:nbox
                val = sqrt(max(variances[offset + box], 0.0))
                @printf(io, "%11.4E%s", val, box == nbox ? "\n" : " ")
            end
        end
    end

    nx, ny, total_slots = size(flxnee)
    open(joinpath(dir, "flxnee.txt"), "w") do io
        for slot in 1:total_slots
            for iy in 1:ny
                for ix in 1:nx
                    @printf(io, "%11.4E%s", flxnee[ix, iy, slot], ix == nx ? "\n" : " ")
                end
            end
        end
    end
end

function needs_hourly_rescale(units::Union{String,Nothing})
    units === nothing && return true
    lowered = lowercase(strip(units))
    isempty(lowered) && return true
    occursin("h-1", lowered) && return false
    occursin("/h", lowered) && return false
    occursin("per hour", lowered) && return false
    occursin("hour", lowered) && return false
    return true
end

function compute_fortran_prior_variances(cfg::Config, domain::Domain, fluxes::PriorFluxes, date_info::DateInfo)
    nbox = domain.nbox
    nbox == 0 && return nothing

    ntstate = max(date_info.ntstate, 1)
    ndt = max(date_info.ndt, 1)
    nd_nee = date_info.nd_nee
    nee_field = fluxes.nee !== nothing ? fluxes.nee : fluxes.nee_global

    if nd_nee <= 0 || nee_field === nothing
        return nothing
    end

    times = nee_field.time
    isempty(times) && return nothing

    statres_int = max(1, Int(floor(cfg.statres + 1e-6)))
    nday = max(date_info.nday, ntstate * statres_int)

    slot_hours = cfg.nstep_nee_reg > 0 ? cfg.nstep_nee_reg : (cfg.nstep_nee > 0 ? cfg.nstep_nee : (nd_nee > 0 ? Int(round(HOURS_PER_DAY / nd_nee)) : 24))
    slot_hours = max(slot_hours, 1)
    slot_step = slot_hours / HOURS_PER_DAY

    total_slots = nd_nee * nday
    nx = domain.nxregrid
    ny = domain.nyregrid

    flxnee = zeros(Float64, nx, ny, total_slots)
    for k in 1:total_slots
        target = date_info.juldatei + (k - 1) * slot_step
        idx = nearest_time_index(times, target)
        flxnee[:, :, k] .= nee_field.data[:, :, idx]
    end

    hour_scale = needs_hourly_rescale(nee_field.units) ? SECONDS_PER_HOUR : 1.0
    scaling = nee_field.scaling_factor / hour_scale
    flxnee .*= scaling

    reg_lat = !isempty(domain.reg_lat) ? domain.reg_lat : [domain.rlly + (iy - 1) * domain.rdy for iy in 1:ny]
    nbox_xy = domain.nbox_xy
    area_box = domain.area_box
    tmp_box = zeros(Float64, nbox)
    reg_nee = zeros(Float64, nx, ny, nd_nee)
    nee_box = zeros(Float64, nbox, ntstate, nd_nee)

    for state in 1:ntstate
        fill!(reg_nee, 0.0)
        for day_idx in 0:statres_int-1
            n = (state - 1) * statres_int + day_idx + 1
            start_slot = (n - 1) * nd_nee + 1
            stop_slot = min(start_slot + nd_nee - 1, size(flxnee, 3))
            if stop_slot >= start_slot
                reg_nee[:, :, 1:(stop_slot - start_slot + 1)] .+= @view(flxnee[:, :, start_slot:stop_slot])
            end
        end
        reg_nee ./= statres_int

        for slot in 1:nd_nee
            fill!(tmp_box, 0.0)
            for jy in 1:ny
                area = grid_area(reg_lat[jy] + 0.5 * domain.rdy, domain.rdy, domain.rdx)
                for ix in 1:nx
                    box = nbox_xy[ix, jy]
                    if box > 0
                        tmp_box[box] += reg_nee[ix, jy, slot] * area
                    end
                end
            end
            for box in 1:nbox
                area = box <= length(area_box) ? area_box[box] : 0.0
                nee_box[box, state, slot] = area > 0 ? tmp_box[box] / area : 0.0
            end
        end
    end

    err_box = zeros(Float64, nbox, ntstate, ndt)
    flxerr = cfg.flxerr
    hour_per_slot = nd_nee > 0 ? 24.0 / nd_nee : 24.0
    bin_width = ndt > 0 ? 24.0 / ndt : 24.0
    hloc_len = length(domain.hloc)
    counts = zeros(Int, ndt)

    for state in 1:ntstate
        for box in 1:nbox
            fill!(counts, 0)
            for bin in 1:ndt
                err_box[box, state, bin] = 0.0
            end
            hloc_offset = box <= hloc_len ? Int(trunc(domain.hloc[box])) : 0
            for slot in 1:nd_nee
                base_hour = (slot - 1) * hour_per_slot
                local_hour = mod(base_hour + hloc_offset, 24.0)
                bin = ndt > 0 ? clamp(Int(floor(local_hour / bin_width)) + 1, 1, ndt) : 1
                err_box[box, state, bin] += nee_box[box, state, slot] * flxerr
                counts[bin] += 1
            end
            for bin in 1:ndt
                cnt = counts[bin]
                if cnt > 0
                    err_box[box, state, bin] /= cnt
                else
                    err_box[box, state, bin] = 0.0
                end
            end
        end
    end

    variances = zeros(Float64, nbox * ndt)
    eps_floor = 1e-12
    for bin in 1:ndt
        mean_vals = sum(err_box[:, :, bin]; dims=2) ./ ntstate
        idx_start = (bin - 1) * nbox + 1
        variances[idx_start:idx_start + nbox - 1] .= max.(abs.(vec(mean_vals)), eps_floor) .^ 2
    end

    debug_dir = get(ENV, "FLEXINVERT_DEBUG_OUTPUT", "")
    if !isempty(debug_dir)
        try
            write_debug_prior_diagnostics(debug_dir, nee_box, err_box, variances, flxnee, nbox, ntstate, nd_nee, ndt)
        catch err
            @warn "Failed to write prior diagnostics" directory=debug_dir err
        end
    end

    return variances
end

function determine_time_steps(cfg::Config)
    date_info = initialize_dates(cfg)
    return max(date_info.ndt, 1)
end

"""
    determine_n_modes(eigenvals::Vector{Float64}, cfg::Config) -> Int

Determine number of modes to retain based on configuration.
"""
function determine_n_modes(eigenvals::Vector{Float64}, cfg::Config)
    total_var = sum(eigenvals)
    total_var <= 0 && return length(eigenvals)

    max_eval = maximum(eigenvals)
    max_eval <= 0 && return 1

    # Number of modes from cumulative-variance threshold (defaults to all modes)
    var_threshold = hasfield(typeof(cfg), :variance_threshold) ? cfg.variance_threshold : 1.0
    n_by_variance = length(eigenvals)
    if var_threshold < 1.0
        cumulative = 0.0
        n_by_variance = 0
        for (i, val) in enumerate(eigenvals)
            cumulative += val
            n_by_variance = i
            if cumulative / total_var >= var_threshold - 1e-6
                break
            end
        end
        n_by_variance = max(n_by_variance, 1)
    end

    # Modes surviving the Fortran-style truncation threshold (relative to max eigenvalue)
    trunc_frac = hasfield(typeof(cfg), :trunc) ? cfg.trunc : 0.0
    n_by_trunc = nothing
    if trunc_frac > 0 && isfinite(trunc_frac)
        cutoff = max_eval * trunc_frac
        n_candidates = count(x -> x >= cutoff, eigenvals)
        if n_candidates > 0
            n_by_trunc = n_candidates
        end
    end

    # Positive eigenvalues safeguard (1e-12 relative to max, mirroring Fortran tolerance)
    positive_cutoff = max(max_eval * 1e-12, eps(Float64))
    n_positive = count(x -> x >= positive_cutoff, eigenvals)
    n_positive = max(n_positive, 1)

    n_modes = n_by_variance
    if n_by_trunc !== nothing
        n_modes = n_by_trunc
    elseif var_threshold >= 1.0
        n_modes = n_positive
    end

    return clamp(n_modes, 1, length(eigenvals))
end

end # module PriorCovariance
